Bogoliubov dynamics of condensate collisions using the positive-P representation 
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We formulate the time-dependent Bogoliubov dynamics of colliding Bose-Einstein condensates in terms of a 
positive-P representation of the Bogoliubov field. We obtain stochastic evolution equations for the field which 
converge to the full Bogoliubov description as the number of realisations grows. The numerical effort grows 
linearly with the size of the computational lattice. We benchmark the efficiency and accuracy of our description 
against Wigner distribution and exact positive-P methods. We consider its regime of applicability, and show 
that it is the most efficient method in the common situation - when the total particle number in the system is 
insufficient for a truncated Wigner treatment. 



CZ2 

wo; 
^— > 



-a 
c 

o 
o 



> 

CO 



in 
o 



X 



I. INTRODUCTION 

The collision of two Bose-Einstein condensates (BECs) - 
if the relative velocity is sufficiently high - leads to the for- 
mation of a halo of scattered atoms. This phenomenon has 
been the object of numerous experimental lfH4l7l and theo- 
retical investigations IIU [lfl[l7l - [34ll . The atoms forming the 
halo could be used for precision measurements ll35ll . inter- 
ferometry |0, l36U39ll . or tests of quantum mechanics iKoll . 
Condensate collisions are also related to such phenomena 
as molecular dissociation ll4ll - l56ll . atomic four-wave mix- 
ing Irj [l4ll57l - l60ll . superradiant scattering ll6lll70ll . atomic 
parametric down conversion ll7ll - l77ll . and impact of a BEC 
on a barrier ||78| - |8T| . 

Recently in a series of experimental studies |Q1 0], a 
quantitative analysis of the supersonic collisions of two 
Bose-Einstein condensates was presented. It was based on 
stochastic Bogoliubov equations for a particle field inter- 
acting via a contact potential. In this manuscript we pro- 
vide the details of that method. It relies on solving a set 
of stochastic equations in a plane wave basis, rather than 
a diagonalization of the Hamiltonian. We have found this 
approach to be more effective as it allows one to study 
large scale multi-mode problems that would not be possible 
with direct diagonalization. This is because in phase-space 
stochastic methods, such as presented in this work, the com- 
putational requirements (memory, time) scale linearly with 
the number of modes or grid points. 

Several stochastic methods have been used with suc- 
cess in the past to study the scattered atoms in these sys- 
tems. They treated the full atom - field system - in 
contrast to a Bogoliubov expansion applied here - using 
the truncated Wigner fl23[ u&l uM |33p and the positive- 
P representations lfl7i |28[ l29l l32i 13311 . However, these 
are not suitable for a majority of current experiments, in- 
cluding the recent metastable Helium condensate collisions 
H S 111 [H [01 [H. The truncated Wigner approach is 
limited to the case when the total number of atoms in the 
system is much larger than the number of necessary modes 
lr2q l82l 18311 . otherwise significant discrepancies ("trunca- 
tion") with full quantum dynamics appear. The positive-P 



approach is complete, but has numerical instabilities that 
make it useful only for short times ll28ll33ll84ll . often shorter 
than the duartion of the collision. 

Instead of a full atom - field approach, a wide class of 
collisions is described accurately by a Bogoliubov descrip- 
tion. This approach is valid while the number of particles 
scattered during the collision is small in comparison with 
the total, a condition satisfied in most of the experiments. 
The time-adaptive refinement, where the condensate wave 
function undergoes mean-field evolution, is sufficient to de- 
scribe most collision experiments. Moreover, contrary to a 
common fallacy, the Bogoliubov formulation takes into ac- 
count the later Bose enhancement and stimulated scattering 
into quasiparticle modes that can occur. 

The drawback of the Bogoluibov method has been that 
accurate description of the real experimental situation typi- 
cally requires a computational grid with 10 6 — 10 7 points. A 
major contributing factor to this large lattice size is the need 
to resolve the supersonic wavelengths in the whole collision 
region. This large lattice renders a direct solution of the 
Bogoliubov-de Gennes evolution equations impossible. To 
avoid the diagonalization, one can introduce a phase space 
distribution for the Bogoliubov field. We call this approach 
stochastic time-adaptive Bogoliubov (STAB). Here we use 
a positive-P representation of the scattered particles, which 
differs from a previous well-known stochastic formulation 
l85ll . which used a Wigner representation. As is demon- 
strated below, the advantage of the present method is a much 
better signal-to-noise ratio in the calculations for the most 
typical regimes of interest. As our positive-P based method 
bases on the broken-symmetry Bogoliubov description, it is 
applicable when the scattered particles are well separated in 
momentum-space from the condensates. This is the case 
for a wide range of supersonic phenomena, which apart 
from the condensate collisions include molecular dissoci- 
ation, superradiant scattering, and parametric down conver- 
sion, as well as supersonic flow past barriers and other im- 
purities JsSH. 

The paper is organized as follows. Section|II]provides the 
Bogoliubov description of a BEC collision. Section [III] in- 
troduces its positive-P representation, and describes the re- 
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suiting stochastic evolution equations used for simulations. 
In Section [TV] we compare the accuracy and efficiency of 
this positive-P Bogoliubov method (P-STAB) with the prior 
trunacted Wigner, positive-P and Wigner Bogoliubov (W- 
STAB) methods for several characteristic BEC collision ex- 
amples. We conclude with SectionM 



The derivation of above equation is standard, and based on 
removing higher-order dependence on 6 and - equiva- 
lent to assuming that the influence of the Bogoliubov field 
on itself is negligible as compared to the impact of the con- 
densate. 

The initial state of the trapped BEC is a solution of the 
stationary GP equation 



II. COLLIDING CONDENSATES - THE BOGOLIUBOV 
DESCRIPTION 

We consider a zero-temperature, single-species bosonic 
gas. As it is dilute, the interatomic interaction can be effec- 
tively reduced to a contact delta potential with strength g. 
In the second quantization, the Hamiltonian reads 

H = Jd 3 x&(x) (-|^V 2 + V(x)^ #(x) 

+ | y<i 3 x* t (x)* t (x)*(x)*(x), 

where m is the atomic mass, g = A-Kh 2 a s /m with a s being 
the s-wave scattering length and V(x) is the external trap- 
ping potential. The field operator ^(x) anihilates an atom at 
position x and satisfies the bosonic commutation relations. 

In order to start the (half-) collision, a superposition of 
two counter-propagating mutually coherent atomic clouds 
is prepared by a Bragg pulse. Simultaneously, the trapping 
potential is turned off. The two fractions start to move apart 
along the z axis with relative speed 2v rcc , twice the atomic 
recoil velocity. We define speed of sound using the den- 
sity at the center of the initial condensate (n ma x), obtain- 
ing c max = \J gn max / m. In the supersonic limit, when 
-max, the gas is no more superfluid and a certain 
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portion of atoms is scattered incoherently out of the BECs, 
forming the halo. The main focus of experiments and the- 
ory are the properties of the atoms in this halo. 

In the time-dependent Bogoliubov approach (we use the 
simpler U(l) symmetry -breaking variety), the field operator 
is split into 



*(x,t) = 0(x,i) + (5(x,i), 



(1) 



where </>(x, t) is the condensate wave function normalized 
to N - the number of particles. Its dynamics is governed by 
the Gross-Pitaevskii (GP) equation 
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d(j)(x, t) 
dt 



_^ y2 
2m 



(x,t). (2) 



The Bogoliubov field operator <5(x, t) describes the "non- 
condensed particles", and obeys the equation 
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with chemical potential /i. The Bragg pulse transforms the 
condensate wave-function into 
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where fco = mv ICC /his the wave-vector associated with the 
recoil velocity. Neglecting quantum depletion, which is tiny 
in most cases, the state of the non-condensed particles is a 
vacuum, denoted by 1 0) . 

A common approach now would be to diagonalize the 
equation (0 using a Bogoliubov transformation, and solve 
the obtained Bogoliubov-de Gennes equations. However, 
for many systems of interest it requires 10 6 — 10 7 points in 
space x, which prohibits such diagonalization. 

Instead, we develop an equivalent stochastic description 
of equation (0 using the positive-P representation. To ob- 
tain the dynamical equations, it is necessary to start from a 
Hamiltonian description. The equation (0, together with its 
conjugate, can be used to trace back the effective Hamilto- 
nian for the Bogoliubov field, 



H cS = J d 3 x^( x ) (-|^V 2 ) 5(x) 
+2g f d 3 x|0(x)| 2 5t( x )5(x) 



(6a) 
(6b) 



:0(x) 2 5 t (x)<5 t (x) + h.c. (6c) 



The line ( I6ab contains the kinetic energy of the noncon- 
densed particles and j6bl the interaction between conden- 
sate and noncondensate particles. Finally, (l6cb governs the 
transfer of atomic pairs from the BEC to the 5 field. 



III. STOCHASTIC TIME-ADAPTIVE BOGOLIUBOV 
(P-STAB) METHOD 

A. Positive-P representation of the Bogoliubov field 

We employ the positive-P representation to expand the 
density matrix for the uncondensed field j(x, t) as a distri- 
bution P over local coherent states at each point x in space, 



p= IP ip,tp A tp,tp V 2 ipV 2 ip 



(7a) 



3 



where the complex fields ip(x) and i/j(x.) are the amplitudes 
of the local off-diagonal coherent state projectors A 



A V,V>J =Q9A x ^(x),V(x)J . (7b) 

X 

Af e ^ ^ x ^ t ( x ) dx |0)(0|e^'^ x )*'^ x)dx . 



with normalisation J\f = e / 'A( x )*'/'( x ) The operator 
|0)(0| projects onto the vacuum state. As the numerical 
computation is made on a grid, the local projectors A x take 
on the form 



^ =e *(x)i'(x)AV 



,«x)*(j(x)-</,(x))AV (?c) 

(7d) 



where A V = Ax ■ Ay ■ Az is the volume per grid point, a = 
?/>(x)V AV, a = ^(x)v 7 AV, and \a) x is a coherent state at 
location x with complex amplitude a. We underline that the 

distribution P ip,ip contains complete information about 
the density matrix p. 

Since it is non-negative and real, it can be regarded as a 
probability distribution of the complex valued fields if>(x) 
and ^(x). It is therefore also equivalent to a large ensemble 
of samples of the fields. Consequently, the state p is re- 
produced by the set of ip(x) and ip(x) when the number of 
samples S tendts to infinity. The assumption that the initial 
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state of 5 is vacuum translates into 

^(x,0)=^(x,0) = 0. 

B. Dynamics 

The quantum evolution of the state 

H e s,p 



(8) 



.,dp 



(9) 



is equivalent to a partial differential equation for P [8S 
IsHUDl, which can be derrived using the operator identities 



(5(x)A = t/>(x)A 



* t (x)A = 
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A<5(x) = 
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(10) 
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These identities are used to convert the quantum operators 
in H c f! inside Eq. (O to partial derivatives. The resulting 
equation is of a Fokker-Planck type and it is well known 
that it can be rendered into a random walk of the samples of 
■0(x) and tp(x) - the Langevin equations - which in the Ito 
representation read 



(14a) 
(14b) 



V 2 + 2 5 |<Xx, t)\ 2 \ t) + g 0(x, i)>(x, t)* + ^ihg 0(x, i)£(x, t) 



, t) + g </>(x, t) V(x, t)* + </>(x, t)£(x, t). 



Here £(x, t) and £(x, £) are delta-correlated, independent, 
real gaussian stochastic noise fields with variances 

(£(x,t)£(x',0) =0 (15) 
(£(x, i)e(x', t')> = i)C(x', *')) = -5 (3) (x - x')<5(i - 0- 

and zero mean. Numerically, £ and £ are usually ap- 
proximated by real gaussian random variables of variance 
l/(AtAV) that are independent at each point at the com- 
putational lattice, and at each time step of length At. 

One important feature of the equations (fl4l is that, simi- 
larly to Eq.©, they are linear in tp and i/j. This way, the 
nonlinear instabilities are absent, together with boundary 
term systematics l92l l93ll and finite-simulation-time issues 
ll84ll that may occur in direct positive-P treatments of the 
full boson field ^. 



C. Observables 

Expectation value of any normal-ordered observable are 
evaluated using the positive-P representation by substituting 
5^ — )• ip* and 6 — > ip and calculating a stochastic average 
JH, denoted as ( ) st , 



( n ^ n *c*)) = & < n n ^) st 

j k j k 

(16) 

Note that since Eq.@ is linear, and the initial state is vac- 
uum, then at all times t 



(0|*(x,t)|0)=0. (17) 



As an example, the one-particle density matrix is given by 

Pl (x,x',t) = (*t( X) ^ (x ' )t )^ 

= 0(x,i)>(x',t) + ( ( 5 t (x,i)<5(x',O> 
+0(x^)*(5(x^> +0(xV)(<5t(x,t)) 

= <f>(x,t)*<f>(x',t) + Wx,f)^(x',i)) lt 

where we used (T7\ on the second-last line. The number of 
non-condensed atoms is 

SN = ( (<5t(x)5(x)) d 3 x = / (V7(x)^(x))) st d 3 x. 



(18) 

For a general observable F, the best estimate of its expec- 
tation value (F) is given by the mean of its corresponding 
estimator f(ip,ip,(f>) 

F = (/> st . 

The uncertainty in this mean is best estimated via the vari- 
ance of a set of subensemble means: We divide the S re- 
alizations into n bins of equal size s (so S = sn), and 
the jth subensemble (j = l,...,n) gives a subensemble 
mean Fj = (f) s t,y Due to the central limit theorem, these 
subensemble means are approximately normal distributed 
(which is not necessarily the case for the estimators from 
individual realizations). As a result, the uncertainty in the 
final mean (F = ^ Y^j=i Fj a ls°) is well estimated by 



'var F 



n- 1 



(19) 



D. Orthogonality and applicability 

It is well known that the U(l) symmetry breaking Bo- 
goliubov method reveals some problems at longer evolution 
times. These are related to an incomplete treatment of the 
the phase spreading of the condensate [r9411 . As the approach 
does not preserve the orthogonality of the non-condensed 
field <5(x) to the condensate mode 19511 . the part of 6 that 
accumulates atop the condensate could just as well be con- 
sidered to still be part of the BEC, and discounted from the 
number of scattered particles. 

For this reason, the results of the above method should be 
treated with caution when modes having significant overlap 
with the condensate are relevant. In practice, such modes 
lie in parts of k-space close to the condensate clouds. For- 
tunately, the bulk of the halo is well separated from the con- 
densates and remains unaffected. 

More generally, supersonicity always leads to orthogo- 
nality between scattered and condensed atoms because the 
condensate mode function contains no plane-wave compo- 
nents above the speed of sound. This allows the use of the 
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FIG. 1. The convergence of observable estimates in the two Bo- 
goliubov methods as the number of trajectories is increased. The 
system is the 4 He collision of (J]. The quantity shown is the to- 
tal number of atoms in the halo at t = 120/xs, well after the end 
of the collision. Narrow k-space regions containing the conden- 
sates were excluded from the atom sum. Black solid line: Wigner 
Bogoliubov calculation (W-STAB), Red dashed: Positive-P Bo- 
goliubov (P-STAB). 

method presented here for collisions of BECs, molecular 
dissociation, superradiant scattering, parametric down con- 
version, or flow past barriers and other impurities. 



IV. RELATIONSHIP WITH COMPARABLE METHODS 

Stochastic evolution equations have been previously de- 
rived for Bogoliubov descriptions of cold atoms systems by 
Sinatra et aJ. l85ll using the Wigner representation. An im- 
mediate question is how the positive-P based method pre- 
sented here compares. We expect that the positive-P method 
will tend to be inherently less "noisy" initially due the lack 
of starting noise which is necessary to represent the vacuum 
in the Wigner treatment. It is also instructive to compare 
performance and accuracy with the two other stochastic 
methods used previously (positive-P and truncated Wigner) 
which treat the whole atom field ^> as one unit without us- 
ing the Bogoliubov approximation. In this section we will 
benchmark these four simulation methods. 



A. Wigner STAB 

Representing the U(l) symmetry breaking description of 
Sec. HI] using the Wigner representation we obtain the fol- 
lowing stochastic description of the field S (x). There is only 
one complex field ip w (x), with the initial vacuum described 
by a random initial condition that places half a virtual parti- 
cle into each mode 



lMx,0) 



£(x,0)+*£(x,0) 



(20) 
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FIG. 2. Slices of the halo density on the plane k z = 0, perpendicu- 
lar to the collision direction, for the 4 He collision of 
right at the end of the collision. Both results are from ensembles 
of 224 realizations, (a): Wigner Bogololiubov, (b): Positive-P Bo- 
goluiubov. 

(The noises £ and £ are as defined by < fT3T > ). The subsequent 
evolution contains no noise and is 



it) 



dip w (^t) 
dt 



2m 



V 2 +2«7|0(x,i)| 2 U w (x,t) 



+gcf){x,t) 2 il) w fx,*)* 



(21) 



Observable calculations differ somewhat because the half- 
particle occupation of the initial modes must be corrected 
for. For example, 



pi(x,x',t) 



(22) 



(x, t)*0(x', f) + (^(x, i)> ro (x', t)) st - -<5(x - x'). 



This is the symmetry-breaking analogue of the more in- 
volved number-conserving description of Sinatra et aJ. l85ll . 
and shares the same noise properties. However the same 
orthogonality caveats (Sec. lIIIDb apply as for the P-STAB 
method derived in this paper. 



B. Full-field methods 

For some parameters, another good alternative is to use 
the truncated Wigner representation to simulate the com- 
plete boson field directly, as was done by Norrie et al. 
j23l HfJ]. This has the advantage of being applicable be- 
yond the undepleted source approximation. However, the 
total number of particles should be significantly larger than 
the number of modes (for correctness ll26l I28II ). This ap- 
proach requires the truncation of some high-order terms in 
the partial differential equation for the resulting phase space 
distribution P, leading to the name "truncated" Wigner rep- 
resentation. Here there is one complex field ?/>n/(x) (no 
separate condensate field 0(x, t)) with the initial state 



V>w(x,0) =0(x,O) 



£(x,0)+z£(x,0) 



(23) 



The subsequent evolution contains no noise and is 

h 2 



dip w (x,t) 

in ; 

dt 



^V 2 + 5 |tMx,;)| 2 ^w(x,t). (24) 



The one-particle density matrix is given by 

Pi(x,x',t) = (^(x,t)*Vw(x / ,t)) s t-^(x-x'). (25) 

Finally, a direct treatment of the full field u sing the 
positive-P representation has been used[ 17, 28, 29[ l32l l33ll . 

Here there are two complex fields Vv( x ) anc ^ Vv( x ) w i tn 
the initial state 



V; P (x,0) = V P (x,0) 

The evolution is 

_dip p (x., t) 



(x,0). 



(26) 



iti- 



dt 



2m 



V- +# p (x,t)> p (x,i) 



(27) 



it) 



dip p (x., t) 
dt 



-— V 2 + #p(x,i)>p(x,i) 



ihgt(x,t)}l> p (x,t). (28) 
The one-particle density matrix is calculated with 

Pi(x,x',t) = $ p (x,t)*^ p ( X ',t)) st . (29) 



C. Efficiency measures 

When considering the halo, the most pertinent observ- 
ables have been the total number of particles, the density 
distribution in k-space, and density correlations between 
specified regions in the halo. Accuracy in the latter two 
kinds of observables hinge on a good signal-to-noise ratio 
of the local density in k-space. The uncertainty of the final 
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estimates is given by ( fT9] l, a function of the ratio between 
variance of the estimator and the number of realizations 
S oc n. So, other things being equal, the computational 
effort required to achieve a set accuracy will scale as that 
variance. Accordingly, in Figs. [3] and |4] (upper panels) we 
will show how the variance of the estimators of halo density 
in k-space compare between methods as a function of time. 
In Figs.Q]and[2]we directly show the noise that is seen with 
the full-field positive-P and Wigner Bogoliubov treatments. 



D. The low and moderate particle number case 

Let us first consider the common case when the total 
number of atoms in the halo is quite low - so low that the 
number of halo atoms per mode is much less than one. Here 
we expect the initial noise in the Wigner methods d20t or 
d23l to be a severe problem, since the initial atom number 
variance there is 1/2 per mode, regardless of how many true 
atoms are present. 

The first plots (Figs. Q] and [2]i are from Wigner and 
positive-P Bogoliubov simulations using the experimental 
parameters of which described the collision of a BEC 
of metastable 4 Hc* atoms. They show the amount of noisy- 
ness in observables after the end of the collision. In this 
case, no bosonic enhacement of the scattering process oc- 
cured, thus the total number of atoms in the halo was quite 
low («1300), while the number of modes was 2.95 x 10 6 . 

The next figure, [3] shows the halo density variance and 
the total number of scattered atoms in the collision of a BEC 
of 150 000 23 Na atoms. This case was considered in several 
previous works 11281,1321, 13311 . Here the halo reached 1.1 x 10 4 
atoms with 1.08 x 10 6 modes). 

We see that the noise in the Wigner calculations is se- 
vere in these cases, as compared to the positive-P methods. 
Although the noise in the P-STAB calculation grows with 
time, it never surpasses the level of the Wigner methods for 
the timescales shown. The variance in both Wigner methods 
is identical. 

The lower panel of Fig.[3]shows the accuracy of the meth- 
ods. Both Bogoliubov methods agree perfectly with each 
other, and with the exact calculation that uses the positive-P 
representation of the full field (for as long as it lasts). The 
truncated Wigner displays a false growth of the number of 
particles in the halo. This is due to known spurious scatter- 
ing by virtual particles when the momentum cutoff is this 
large, as described in i28l l33i l82ll . For this simulations, the 
number of spatial modes is much larger than the number of 
true particles (150 000). 



E. The high particle number case 

A different situation is presented in (Fig. [3}, where we 
used the parameters from 112311 . where 6 x 10 6 atoms of 23 Na 
participated in the collision. There were 3.14 x 10 6 spatial 
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FIG. 3. 23 Na BEC collision as in HHHIH. N = L5 x 1C|5 
Upper panel: Variances of local atom density estimators in the 
slice at k z = obtained for various methods, as for use in (1191 
- see text. An average value over all k x and k y locations in the 
slice is shown. Lower panel: Number of scattered atoms in the 
halo. Solid red: Positive-P Bogoliubov simulation as described 
in this paper; Blue circles: Wigner Bogoliubov simulation; Dot- 
dash green: truncated Wigner simulation; Dashed black: positive- 
P simulation of the full field. 

modes. As final depletion of the condensate is large (about 
40%), and the Bogoliubov calculation must was stopped at 
t w 280/is, when depletion was 10%. Indeed, in the lower 
panel of Fig [3] one sees a difference beginning to appear 
between the two simulations at this time. In comparison, 
significant dynamics lasts until w 1000/is (not shown). 

The noise performance of the P-STAB method is superior 
here only for t < 300/xs. However, this still matches the 
entire period when the Bogoliubov description is accurate. 



V. CONCLUSIONS 

We have developed the above positive-P Bogoliubov 
stochastic simulation method for use with cold atom gases 
and benchmarked it with existing approaches. As with other 
phase-spae methods, it lends itself to simulation of quite 
general systems, as the calculation is carried out on a sim- 
ple rectangular grid in x / k space, and individual realiza- 
tions are run independently of each other. The computa- 
tional complexity involved scales linarly with the size of the 
computational lattice used, allowing for up to w 10 7 points 
in the lattice on a common workstation. 

The method is applicable for a wide range of supersonic 
phenomena, its main limitations being (1) that the bulk of 
scattered atoms are well separated from the condensates in 
momentum space, and (2) that the depletion of the original 
condensates can be neglected. The condensate wave func- 
tion is, however, free to evolve in time. We note particularly 
that the method handles both spontaneous and stimulated 
scattering. 
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FIG. 4. Na BEC collision with the same parameters as used in 
l23ll . N — 6 x 10 6 Upper panel: Variances of local atom density 
estimators in the slice at k z — obtained for various methods, as 
for use in dl9t - see text. An average value over all k x and k y 
locations in the slice is shown. Lower panel: Number of scattered 
atoms in the halo. Solid red: Positive-P Bogoliubov simulation as 
described in this paper; Blue circles: Wigner Bogoliubov simula- 
tion; Dot-dash green: truncated Wigner simulation; Dashed black: 
positive-P simulation of the full field. The Bogoliubov simulations 
were stopped when the depletion reached 10%. 

The positive-P Bogoliubov method is superior in effi- 
ciency to the Wigner representation in almost all cases that 
we have seen where a U(l) symmetry breaking Bogoli- 
ubov method can still be applied. However, one can imag- 
ine some long time situations where the Wigner simula- 



tion wins since, other things being equal, the variance in 
the positive-P approach grows approximately linearly with 
time, while the variance in the Wigner method stays ap- 
proximately constant around its initial, large, value (These 
trends are seen in the top panel of Fig. [3] ). For situa- 
tions where the overlap between the scattered and conden- 
sate field is non-negligible the number-conserving Wigner 
method l85ll can be used instead. For situations with large 
condensate depletion, there remain the truncated Wigner or 
positive-P treatments of the full boson field. 

A more robust positive-P formulation that explicitly 
imposes orthogonality between condensate and quasi- 
particle modes as in the number-conserving Bogoliubov 
treatment ll95ll is under development and will be presented 
in a forthcoming work. 
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